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On the transient models of the VITAS code: applications in the C5G7-TD pin-resolved benchmark 


problem 
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This article describes the transient models of the neutronics code VITAS that are used for solving time- 
dependent, pin-resolved neutron transport equations using finite elements. VITAS uses the stiffness confine- 
ment method (SCM) for the temporal discretization in order to transform the transient equation to a transient 
eigenvalue problem (TEVP). To solve the pin-resolved TEVP, VITAS incorporates the heterogeneous varia- 
tional nodal method (VNM). The spatial flux is approximated at each Cartesian node using finite elements in 
the z — y plane and orthogonal polynomials along the z axis. The angular discretization employs an even-parity 
integral approach in the nodes and spherical harmonic expansions at the interfaces. To further lower the compu- 
tational cost, the predictor-corrector quasi-static SCM (PCQ-SCM) is developed. Within the VNM framework, 
computational models for the adjoint neutron flux and kinetics parameters are presented. The direct-SCM and 
PCQ-SCM are implemented in VITAS and validated using the 2D and 3D exercises of the OECD/NEA C5G7- 
TD benchmark. In 2D and 3D problems, the discrepancy between the direct-SCM solver's results and those 
reported by MPACT and PANDAS-MOC is less than 0.97% and 1.57%, respectively. In addition, numerical 
studies comparing the PCQ-SCM solver to the direct-SCM solver demonstrate that the PCQ-SCM allows for 
a substantially larger time-step size, hence reducing the computational cost by a speedup factor of 10 to 100 


without compromising numerical accuracy. 
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I. INTRODUCTION 


Neutronic calculations play a critical role in analyzing 
different kinds of nuclear systems[1-3]. Without spa- 
tial homogenization and the application of lower order 
approximations[4], high-fidelity neutron transport methods 
can accurately predict steady-state and transient behaviors 
of neutron flux and thermal power during operation or se- 
vere accident in nuclear reactors. Solving the transient pin- 
resolved equation is particularly challenging since the cal- 
culation is computationally expensive, and it is difficult to 
reach a balance between acceptable numerical accuracy and 
computational costs. The method of characteristics (MOC) 
and the finite element method (FEM) have shown promise 
in this area for high-fidelity calculations. In recent years, 
a number of high-fidelity neutronics codes have been de- 
veloped. The MOC-based codes MPACT[4], PROTEUS- 
MOC[5], PANDAS-MOC[6], NECP-X[7, 8] and the FEM- 
based code Rattlesnake[9] are typical. 

In addition, the variational nodal method (VNM) has 
also shown tremendous promise in solving high-fidelity 
problems[10]. The VNM is constructed from the func- 
tional form of the second-order even-parity transport equa- 
tion, where odd-parity Lagrange multipliers are applied to 
enforce the nodal balance[11, 12]. It has been demon- 
strated that VNM is compatible with a variety of discretiza- 
tion techniques and node geometries. Due to its adaptability 
and accuracy, the VNM has been successfully implemented 
in neutronics codes including VARIANT[13], VIOLET[14], 
VITAS[15], etc. In particular, VITAS is a multi-purpose neu- 
tron transport solver that has the high-fidelity modeling ca- 
pability (namely heterogeneous VNM). Therein, the spatial 
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flux distribution is approximated with quadratic finite ele- 
ments and orthogonal polynomials. Within the node, the 
integral approach for angular discretization is implemented, 
while spherical harmonics functions are used at the interfaces. 
These discretization approaches enable VITAS to accomplish 
refinements in both space and angle. 


For the performance of high-fidelity calculations in tran- 
sient cases, efficient and accurate temporal discretization 
schemes are essential[16]. The quasi-static approach is the 
most often adopted temporal discretization scheme among 
high-fidelity methods. It consists of an improved quasi- 
static method and a predictor-corrector quasi-static (PCQ) 
approach[17, 18]. In contrast to direct methods, the computa- 
tional efficiency of the PCQ method is accomplished by cou- 
pling the transport equation and exact point kinetics equation 
(EPKE) based on the observation that the flux shape varies 
more slowly with time than the flux amplitude. The shape 
function is calculated by transport calculation with coarse 
time-steps, but the amplitude function is determined by EPKE 
calculation with fine time-steps. Consequently, the high- 
fidelity transport codes utilizing the PCQ approach can pro- 
vide time-dependent solutions at an acceptable computational 
cost. 


In the PCQ method, the transient problem is typically trans- 
formed into a transient fixed source problem (TFSP) based on 
the backward difference method (BDM), and is then solved 
by TFSP solvers. In the absence of an efficient acceleration 
method, however, the TFSP will converge slowly[5]. The 
stiffness confinement technique (SCM) is another method that 
can be utilized for treating the temporal dependence. As op- 
posed to the BDM, the SCM transforms the transient problem 
into a transient eigenvalue problem (TEVP). The method was 
originally proposed to alleviate the stiffness in the point kinet- 
ics equation by incorporating the dynamic frequencies[19], 
and subsequently extended to be used in time-dependent dif- 
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fusion and transport problems[20—22]. 

Recent studies on the use of the VNM to time-dependent 
problems have proved its promise. For instance, the VNM 
with the SCM has been examined in whole-core transport 
calculations with spatial homogenization[23]. It was proven 
that the VNM with SCM can leverage a coarse time-step to 
drastically minimize computing cost. To obtain an accurate 
flux amplitude, However, the traditional SCM (denoted as 
the direct-SCM) still requires finer time-steps than the PCQ 
method. In high-fidelity transient problems, it may impose 
large computing overhead. To extend the accurate and effi- 
cient application of VNM and SCM to high-fidelity transient 
situations, it is necessary to develop an improved VNM-SCM 
framework that incorporates the advantages of other temporal 
discretization techniques, such as the PCQ method. 

In this work, the direct-SCM is employed in conjunction 
with the VITAS code to solve time-dependent high-fidelity 
transport problems. The PCQ-SCM is proposed and imple- 
mented in VITAS to reduce the computational cost, which in- 
corporates the fundamental idea of PCQ, but obtains the flux 
shape via solving TEVP, which is transformed by SCM. The 
purpose of the study is to demonstrate the capability of the 
heterogeneous VNM for solving transient pin-resolved prob- 
lems. The remainder of the sections are structured as fol- 
lows: methods for the direct-SCM, heterogeneous VNM and 
PCQ-SCM, are introduced and derived in Section II. Section 
III describes the OECD/NEA C5G7-TD benchmark and its 
associated 2D and 3D exercises. The results obtained using 
direct-SCM and PCQ-SCM are compared in detail in Section 
IV. Moreover, comparisons to other codes are given. Section 
V presents the conclusions and summary. 


Il. THEORY 


In the SCM, the time-depedent transport equation is first 
transformed to a TEVP by frequency transformation. Then, 
by solving TEVP using the steady solver iteratively, the nec- 
essary quantities for updating frequencies are obtained. In 
VITAS, the steady solver based on the heterogeneous solver 
is applied to solving the TEVP. In the following sections, the 
frequency-transformed equation is derived which yields the 
form of TEVP. Then, based on the heterogeneous VNM, the 
nodal functional of the TEVP is presented and discretized to 
obtain response matrix equations. Next, the flowchart of the 
essential iteration for updating frequencies in the direct-SCM 
is presented. Last, the PCQ-SCM is formulated based on the 
direct-SCM and the EPKE. 


A. Frequency-transformed equation 


The dynamic frequency of the scalar flux is defined as: 


oy (r,t) (1) 


ge 


we (r,t) = 
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in which the scalar flux ¢, (r,t) of group g is given by: 


s.t) = [ d, (r,t) Q) 
where v, (r, Q, t) is the angular flux at position r in direction 
€). Note that to remove factors of 7 from these equations, dQ 
is normalized such that f dQ = 1. Generally, it is supposed 
that the dynamic frequency of the angular flux is isotropic 
such that: 


1 o 
= a a Q,t 3 
Wg (r, ) m (r, , t) ais (r, , ) ( ) 
Then, the time-variant angular flux can be expressed as: 
4 (0,0 9,00) eet!) A 


The dynamic frequency can be further decomposed as 
wg(r,t) = ws,g(r,t) + wr(t)[21]. The flux amplitude fre- 
quency wrp (t) represents a global quantity and is dependent 
only on time; the flux shape frequency ws,g(r,t) is depen- 
dent on the space, time and energy. With this approximation, 
the angular flux can be decomposed using the amplitude term 
and the shape term as: 


Wg (r, 2, t) = P(t) Wg (r,Q, t) (5) 
in which P (t) is the flux amplitude, and 4), (r, Q, to) is the 
flux shape. The time-variant flux amplitude and flux shape 
are evaluated as: 


P(t)= P(to)ede dt' w(t’) (6) 


Wg (r,Q,t) = Wg (r, Q, to) alts dt ws, (r) (7) 


To make the decomposition unique[21], a physics-based 
constraint on the shape frequency is introduced: 


G ^ 
ES Ezo(rst) f dG, (r,Q,t)=P(to) (8) 
g—1 


where, v is the number of released neutrons per fission reac- 
tion, X g(r, t) is the macroscopic fission cross-section. The 
constraint guarantees that the shape frequency affects only the 
flux shape, not the amplitude. Similarly, the precursor con- 
centration frequency of delayed group i is defined as: 


zo E 
Ci(r, t) Ot 


Hi (r, t) Ci (r, t) (9) 
in which C;(r,t) is the precursor concentration at position 
T, Which is supposed to be isotropic commonly. Introducing 
Eqs. 3 and 9 to the time-dependent transport equation[23] 
with isotropic-scattering, the frequency-transformed neutron 


transport equations are obtained: 
WT (t) + wsg (T, t) 
v, (r) 
5i deg (r, t)vs (r, Q, t) 
= XL t)dg(r, t) B G(r, t) 


Vs(r, Q, t) 4- Q- Vus (r, Qt) 
(10) 
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isa Where X5 (r, t) and X,,(r,t) are the macroscopic total and 
15 Scattering cross-sections, respectively. The group source 
156 qg(r, t) is the contributions from scattering, prompt and de- 


157 layed neutrons: 
t) = 2 Legg! (Tt) bg 


g #9 


(r,t) 


xg (r) [1 — 8 (r)] 
158 + mit T 11 
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x SD r,t) by! (r,t) 

E 
s in which X, (r,t) is the macroscopic scattering cross- 
so Section, Xg (r) and xi, (r) are fission spectra of prompt and 
e: delayed neutrons, 6; (r) and 8 (r) are fractions of delayed 
se neutrons, and A; (r) is the delayed constant of the precursor. 
1683 The frequency-transformed precursor equations are given as: 


pi(r,t)C;(r i) Davee r,t)ó, (r,t) 


= ^i (r T 


Introducing the dynamic eigenvalue kp in Eq. 10 and mov- 
) ee ing the frequency term to the right-hand-side, the equation 
167 can be transformed into an eigenvalue problem (EVP) as: 


Q- Vyar, Q, t) + Err, t)r, Q, t) 


(12) 
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168 j (13) 
= Xsg(r.t)og(r.t) +0,(r,t) 
t) = M Dogg’ (06, (r,t) 
9*9 
WT (t) T Wsg (r, t) 
170 A (r) g(r, t) (14) 
+ kx, (r,t) 3 vE ur (r, t)by (r,t) 


g' 
171 where X; (r,t) is the dynamic fission spectrum, which is de- 
172 fined as: 


I 
Xig (r) Bi (r) Ai (r) 

e-5e Y ET 

173 (15) 
174 After being transformed to an EVP problem, Eq. 13 can be 
175 solved with an existing neutron transport solver. The dynamic 
176 frequencies resulting in the dominant eigenvalue equal to 1 
177 are the solutions to Eq. 13. Thus, the dynamic frequencies are 
178 Solved iteratively using power iteration and the secant method 
179 to non-linear equations. This non-linear iteration algorithm 
180 to frequency-transformed equations is named as the k — w 
181 Iteration[23], which is illustrated in Fig. 1. 


Xg (T, t) = Xg 


182 B. Variational nodal formulation 


183 In this section, the discretized functional of EVP is derived 
184 for the VNM. Last, the corresponding response matrix equa- 
185 tions are presented in a multigroup framework. 


186 1. The discretized functional 


In the VNM, the even- and odd-parity flux yt and v^ are 
188 defined, and the corresponding functional in a node v can be 
written with explicit separation of the radial and axial inter- 
faces as[11]: 


ct +)2 1 pt2 
nisal» [a In teas] 
—Lsb* — 20q 


42 / de a db J dnp - (jw 
Ti; yt de 
+2 f aa [ao ( +n,- Quit |z- ) 
191 I 


192 For brevity, r, Q and t in the unknowns are suppressed. In the 
local coordinate, the node is defined in —Az/2 € z € Az/2, 
the planar area is A = AzAy, n, is the outward normal 
to the lateral interfaces extending over the periphery I’, n,+ 
and n,- are the outward normal to the top and bottom axial 
interfaces, respectively. The functional in the whole problem 
domain is the superposition of the functional in each nodal 
volume as F [V * Y7] = $5, Fp [v * , v]. 

Within the node, the spatial distribution of the even-parity 
flux is approximated by: 
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203 The scalar flux can be expanded as: 
a. O(n, t) = f3 (2) 8g" (x,y) 6 (0) (18) 


where the time-dependent scalar flux moments are defined by 

t) = f dw (Q,t). f, (z) is a vector of orthogonal poly- 
nomials defined in the node v. The polynomials are governed 
208 by: 


Jer (z) f- (2) = AzI (19) 


where J is an identity matrix. g (x, y) is a vector of continu- 
ous finite-element basis functions. The triangular and quadri- 
lateral iso-parametric quadratic finite elements (FE)[11] are 
employed to map curved interfaces between materials. Fig. 2 
shows the FE mesh with 32-elements used to describe a fuel 
2:6 pin cell in this paper. 

The odd-parity flux on the axial interface is approximated 


V^ (r, Q,t) ~ fz (EA2/2) yz (Q) & h^ (x,y) x. (t) 


(52/2) y7, (Q) @AT (x,y) x., (t) 
(20) 
in which h (zx, y) is a vector of piecewise constants. y? (Q) 
and y7, (Q) are odd-parity spherical harmonic vectors com- 
222 prised of both sine and cosine functions[24]. The vectors 
y1 (Q) are low-order P, approximations, while y7, (Q) and 
Vso (Q) contain higher-order terms from n + 2 to some larger 
number JV. In this paper, we refer to these as approximations. 
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Update frequency | 
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Calculate precursor 
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Fig. 1. Flowchart of the k — w iteration 


Fig. 2. A finite element mesh for a fuel pin. The orange and blue col- 
ors represent the fuel region and the moderator region, respectively 


226 X, (t) and Xz (t) are the corresponding low- and high-order 
227 moments of the odd-parity flux on the axial interfaces. 


28 The odd-parity flux on the lateral interfaces is approxi- 
22 mated by: 


T 


y? (0) & FY (C) x, (t) 
+ y2, (Q) & f2, (C) x, (t) 
y-uy C=y,z 


V^ (r, 2,1) = 
(21) 


230 


2s Where y = 1,2,3,4 corresponds to the y = z^,z*,y ,y* 
232 lateral interfaces. In VITAS, P. (C) and n. (C) are chosen 
23 as fourth-order polynomials. Similarly, yZ (Q) and y7, (Q) 
234 are low- and high-order P,, approximations, and x, (t) and 
295 Xyo (t) are the corresponding low- and high-order moments 
23 Of the odd-parity flux on the lateral interfaces. 

Inserting the trial functions of Eqs. 17, 18, 20 and 21 yields 


237 


238 the discretized functional: 


F [px x] = f 0 (0.04(, 090.1) 
— $7 (t) I; & F, (t) 6 (t) - 26 (t)! q (t) 
+25_ J dT (Q, t) EJ (Q)x., (t) 
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241 


= I drf.(z)@g(a,y)q (rt) | Q3 


22 in which A and E are coefficients matrices. The matrices 
243 containing integrals over the spatial trial functions are eval- 
244 uated numerically using standard FE techniques. The defini- 
245 tions of these matrices are detailed in the previously published 
24 research[11, 12]. Requiring the discretized functional Eq. 22 
247 to be stationary with respect to variation in w(Q, t) yields: 


A(Q, t)p(Q, t) _ I, e Fo (t) 
= q (t) - E()x; (t) - E«(9)x, (t) 
248 (24) 


249 Note that the terms in Eq. 24 is re-grouped by low-order (with 
250 Subscript l) and high-order (with subscript o) terms. 


2. Nodal response matrix equations 


22 The odd-parity moments x; and x, are defined to be con- 
253 tinuous across interfaces. Also, requiring Eq. 22 to be sta- 
264 tionary with respect to variations in x; and X, yields that the 
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even-parity flux is continuous across the interfaces: 


e - A^ f QET (v (0.0 (25) 


With integration over the angle and matrix operations[11] 
to Eq. 24, the scalar flux can be determined in terms of the 
response matrix equations given as: 


é(tü-V(t«q()-c([;i* (0-3 OH] 20 


The coefficients matrices can be also found in [11, 12]. 3* 
and j correspond to the interfacial expansion moments of 
outgoing and incoming neutron currents, respectively. They 
are defined as: 


j+ O= 1&0 ix qm 


The outgoing currents moments j * satisfy: 


j'(0- B(tq (t - R(t)3 (t Q8) 


To sum up, the flowchart of the iteration in the VNM is 
illustrated in Fig. 3. 


C. Direct-SCM in VNM 


The application of high-order iso-parametric FEs incurs 
significant computational costs and the memory usage. To 
alleviate this burden, the element-wise flux values are stored 
in the code rather than vertex-wise values generally adopted 
in the FEM. The element-wise average scalar flux is defined 
as: 


E(t) =E | dedyh Gs.) g" (2,4) bolt 29) 


where $4 (t) is the segment of (t) corresponding to the com- 
ponents of the vector of orthogonal polynomials f , (z) of or- 
der 0, and & is a diagonal matrix composed of the FE areas. 
h(a, y) is a vector of piecewise functions, which are equal 
to one on the domain of the corresponding FE and zero else- 
where. By averaging over each FE, the number of stored val- 
ues is significantly reduced. Correspondingly, the flux shape 
frequency and precursor concentration are also defined as the 
element-wise quantity. For brevity, the subscription v and k 
denote the average quantity in element k of node v in Sec- 
tion II C. For example, ¢,,,9 (t) is the element-wise average 
scalar flux of group g in element k of node v. 

The calculation procedure of the direct-SCM within the 
framework of the VNM is given in Fig. 4. In the k — w 
iteration, the secant method[21] can be applied to update the 
amplitude frequency. The update formula of the amplitude 
frequency based on the secant method is shown as: 
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moments, out-currents 
and fission source 


Calculate response 
matrices 
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sources 


Update out-currents of 
group g 
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Update flux moments 
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Update flux moments 
of group G 


Flux moments 
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Fig. 3. Flowchart of the iteration in the VNM 
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297 (30) 28 where m denotes the iteration step of the k — w iteration. The 
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normalized element-wise average scalar flux Dux; (tn) are 
obtained by: 


Óvk.g (tn) ex Óvk.g (tn) 
P (t4 1) 


x - 
X, Er Xa VYuksfo (tn) AvkGo,g (tn) 

- (31) 
in which A, is the element's area, dx,9 (tn) is the element- 
wise average scalar flux obtained from the EVP solver di- 
rectly before the normalization. This normalization enforces 
the total power contributed by the normalized flux to equal 
the power of the previous time-step. Therefore, the shape 
frequency obtained from the normalized flux is independent 
of the flux amplitude. Additionally, it is supposed that the 
shape frequency is homogeneous within the element. Thus, 
based on the isotropic and homogeneous approximation, the 
element-wise shape frequency is calculated as: 


«| 


where At, = tn — t, .1, and vj, s,; (tn) is the element-wise 
shape frequency within [tn—1, tn]. 

It is supposed that the amplitude frequency varies linearly 
within [t4 1, tn]: 


1 
Atn 


dvk,g (tn) 
Qvk,g = 1) 


Uvk,S,g (tn) = 


| (32) 


cT (t) za QT (ts 1) (t 
At 

Then, according to Eq. 5, the actual element-wise average 

scalar flux is: 


WT (t) = wr(tn—1) + tn—1) (33) 


wr (in) tor (tn—1) 
duk,g (tn) = Qvk,g (tn-1) e 2 Anio Sg (tn)Ata 


" erp (n) wm (t4 
E "m (tn) B T( irent 2 At, 

(34) 
The first-order analytical precursor integration method is im- 
plemented to calculate the precursor concentration, where the 
fission source is supposed to vary linearly across the time- 
step [tn—1, tn]. The analytical solution to the element -wise 
precursor concentration is then expressed as: 


Crk i (tn) = Ck i (tn—1) eg Ark itn 


tn 35 
4 fapa] Qui (t)e^v*tgt (35) 
tn-1 


in which the element-wise fission source Q,,x (t) is expressed 
as: 


4 Qur(tn) = Qui (t5 1) 


Quk (t) = Qurltn-1) At (t tn—1) 
n (36) 
Qui (tn) = KUvk,fg (tn) Qvk,g (tn) (37) 


The update formula of the element-wise precursor frequency 
is given as: 


Pe —xX a Cyk i(f4) #0 
vk iltn) = * Ouk (tn) vk ; 38 
Dwk,i (tn) | Ü Cnn t) a (38) 
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D. PCQ-SCM in VNM 


To reduce the computational cost without a significant loss 
of the accuracy, the PCQ-SCM is proposed. The main idea of 
the PCQ-SCM is to firstly solve TEVP with a coarse time-step 
to obtain the predicted neutron flux. Then, by determining ki- 
netics parameters (KPs) with the predicted flux, the EPKE is 
solved with a fine time-step to obtain the amplitude function. 
The time-steps specifications are illustrated in Fig. 5. The 
error of the neutron flux resulting from the coarse time-step 
is then corrected with the amplitude function. 


Coarse time-step 


a 
[E lad tn tatı 


Fine time-step 


Fig. 5. Illustration of the coarse and fine time-steps in the PCQ-SCM 


In the following sections, the PCQ-SCM is derived from 
the solution to the adjoint equation with VNM. Then, the ap- 
proach for evaluating the reactivity and kinetics parameters 
(KPs) with the VNM-based adjoint flux is discussed. 


1. Adjoint solution in VNM 


The solution to the adjoint steady transport equation, i.e., 


s» the adjoint neutron flux, is adopted as a weighting func- 


tion in the perturbation problems. For example, the exact 
point kinetics equation (EPKE) is derived from the first-order 
perturbation theory with the adjoint flux as the weighting 
function[25, 26]. To implement the PCQ method in the tran- 
sient calculation, the VNM based adjoint flux should be used 
as the weighting function. 

The VNM transforms the original transport equation to 
a second-order even-parity equation. Since the leakage 
operator in the second-order even-parity equation is self- 
adjoint[27], the adjoint calculation is simplified significantly. 
The forward and adjoint response matrix equations differ only 
in the source term. As a result, forward and adjoint calcu- 
lations can be performed using the same within-group algo- 
rithms and response matrices shown in Fig. 3. The adjoint 
source term can be constructed by switching the group num- 
ber of the scattering and fission cross-sections. 

In the following derivations, an asteroid denotes the so- 
lution to the adjoint equation. The adjoint response matrix 
equations at the initial condition are given as: 


p“ (to) = V (to) a* (to) —C(to) [3^ (to) — j= (to)] (39) 


3^ (to) = B(to)q' (to) + R (to) j (to) — (40) 


where the response matrices are identical to those in Eqs. 26 


377 and 28. The vector of the initial adjoint source q* (to) is de- 
378 fined as: 


ii a” (to) = J drf,(z)eg(s,y) (r,t) (41) 
(r io) = b» Pusat T ,to)ó y (r, to) 
da g #9 (42) 
+ ko vX go(m, to) X xz (m, to) 6%, (1, to) 
F 
I 
w G(r, to) = xg(r) [1 — 8(7)) + M xis (vr) Bi(r) (43) 
jl 


384 Note that the eigenvalue kef in the adjoint equation should be 
35 identical to that in the forward equation. 


2. Exact point kinetics equation 


The EPKE is derived from the transient transport 
equation[28]. Similarly, the neutron flux is factorized into 
the shape function and the amplitude function as: 


Wg (r, Q, t) =T (t) Pg (r, Q, t) 
where T (t) is the amplitude function, and v, (r, Q, t) is the 


shape function. Similarly, to make the decomposition unique, 
the constraint is introduced as: 


G 
* I T = 
EET (r, Q, to) POL (r,Q,t) 2 T 


(44) 


) 390 


P 
2 


(to) 


(45) 
where 4/7 (r, Q, to) is the adjoint angular flux at the steady- 
state. It indicates the integral of the shape function with the 
initial adjoint flux is held constant over time. By inserting 
Eq. 44 into Eqs. 10 and 12, and multiplying the adjoint flux 
on both sides, then integrating over all variable, the EPKE is 
obtained as follows: 


= I 

401 ST (t) = P Tu (t) T (t) + >, Ai (t) Gi (t) (46) 
d B; (t) XY 

403 qj (t) = A(t) T (t) — Xi (t) G (t) (47) 


4oa Where p (t) is the reactivity, 6 (t) and P; (t) are effective frac- 
os tions of delayed neutrons, A; (t) is the effective decay con- 
stant, A (t) is the prompt neutron lifetime, and C; (t) is the 
47 reduced precursor concentration. The detailed formulation of 
48 EPKE is available in [17]. Eqs. 46 and 47 can be solved by 
o» ordinary differential equation solvers. In VITAS, the ordi- 
aio nary differential equation is developed based on the SCM for 
41 EPKE[19]. 
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As shown in Section II B, the even-parity integral method 
implemented in VITAS only leads to a solution of the scalar 
flux. Reconstructing the angular flux requires the extra stor- 
age of temporary matrices about the discretized operators[11] 
and extra algebraic operations of these matrices to the scalar 
flux. Besides, storing the angular flux also increases memory 
usage significantly. Therefore, the scalar flux is employed as 
the weighting function to calculate KPs in VITAS. Note that 
the reactivity is not included in the KPs in this section. The 
details of KPs calculated with the scalar flux @ and ¢t are 
listed as follows: 


419 
420 
421 


422 


= Xl) p (r, to) 
ro -Y fal xed, fol (r : t) à, (r,t) 
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= Rr at Ayk Pik, o (to) 
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428 (48) 
€ fg. (r, to) Bi (r) | 
8 zoe | dag! VM gg (7,0) oy (r,t) 
B, (t) = FO 
= - AVRO ik 9 (to) Buk, a 
OM- 122, Sm X$ ULL sd (D) dus a 
7 F (t) 
426 (49) 
Salo ET) _ T(t) 
- ee FO pg 9 
a - ;o (r) A; (r) C; (r 
Xe (e) = Desert arb (sto) Xig (r) i (r) Ci (r,t) 


Yea S ards (r, to) xig (r) Ci (r,t) 
Xa De ee Ave Pig (fo) Xvk ig o Cuna ©) 


Xa Lev Lk Avda, (to) Xi ig Ci e 
The related research shows that the KPs calculated with the 
scalar flux agree well with those calculated with the angular 
432 flux[5]. However, using the scalar flux as the weighting func- 
433 tion can result in significant deviations in the reactivity[5, 29]. 
434 For example, in the moderation density reduction problems of 
the C5G7-TD benchmark, the scalar flux weighting function 
underestimates the reactivity by up to 15.6%[5]. The differ- 
ent reactivities resulting from the scalar flux are mainly due to 
that the leakage operator €? - V is highly anisotropic. There- 
439 fore, adopting the isotropic approximation to evaluate the per- 
4 turbation of the leakage operator can cause significant devia- 
44 tion to the reactivity. To obtain an accurate reactivity for the 
4:2 EPKE calculation, the reactivity is evaluated[30] as: 


keg (t) - 1 
keg (t) 


asa Where keg (t) is the dynamic eigenvalue, which is defined as 
455 the eigenvalue to Eq. 13 with setting frequencies u and w to 


435 
436 
437 


438 


p(t) (52) 


ass Zeros. Since keg (t) is obtained by solving transport EVP, Eq. 
447 52 is equivalent to calculating the reactivity with the angular 
4 flux. It eliminates the necessity to reconstructing and stor- 
4 ing the angular flux. Additionally, solving for the dynamic 
eigenvalue does not increase computational cost due to that 
the SCM requires solving EVPs iteratively shown in Fig. 1. 
Denoting reactivity evaluation using Eq. 52 as Method 1 and 
453 using the weight of scalar flux as Method 2, the comparison 
454 Of two methods is shown in Fig. 6. The results indicate that 
455 the reactivity evaluation using scalar flux underestimates the 
4 reactivity by an error of -15.0% and -10.0% in TD3-4 and 
47 TD5-1 of C5G7-TD problem respectively, which is consis- 
ass tent with the results reported by [5]. 
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Fig. 6. Reactivity evaluation using different methods in C5G7-TD 
problems 


3. Overall PCQ-SCM scheme 


459 


40 Inthe PCQ-SCM solver of VITAS, the initial adjoint scalar 
ai flux $7 (r, to) and forward scalar flux @, (7, to) are obtained 
s2 by solving initial EVPs. The initial amplitude function T (to) 
45 is determined with P (to). Then, in each coarse time-step 
464 [t4 1, tn], the TEVP and the EPKE are solved in the follow- 


465 ing steps: 


sæ — l. Solve the predicted scalar flux ¢, 4 (T, tn), flux ampli- 


467 tude P, (tn), and precursor concentrations C, ;(r, t) 


using the TEVP in Eq. 13 with k — w iteration at 
tn. In the first iteration-step, the frequencies jz and w 
are set to zero. Then, the eigenvalue obtained in the 
first iteration-step is the dynamic eigenvalue keg (tn) 
required in Eq. 52. 


468 
469 
470 
471 


472 


a3 — 2. Determine the reactivity p (tn) with Eq. 52. 


474 3 . 


475 


Calculate the predicted shape function dbp, g (T, tr) with 
the constraint in Eq. 45 as: 


dp,g (r, tn) = Pp,g (r, tn) 
T T (to) 
G * 
2c Jy dr f do; (r, to) aay org (7, tn) 
(53) 
4. Calculate the KPs at t, with Eqs. 48-51, and solve 
EPKE to obtain the amplitude function T (tn) at the 


fine time-step using the interpolated KPs from the val- 
ues at tn—1 and tn. 


476 


477 
478 
479 


480 


ai 5. Determine the corrected flux amplitude P.(t,) and 


482 scalar flux $.,, (T, tn) as follows: 


T (tn) P (t 

= P.) = TIRES P 
P, (tn) 

"t e.g (T tn) P, (tn) P” glo) (55) 


P; (tn) Das (r, tn) 


6. Re-calculate the precursor concentrations with the cor- 
rected scalar flux by Eq. 35. 


486 


487 


488 Then, iterate through Step 1 to Step 6 in the next coarse time- 
489 Step until the end of the transient event. The overall flowchart 
490 Of the PCQ-SCM is shown in Fig. 7. 

Additionally, in the PCQ-SCM solver, the amplitude fre- 
42 quency wry is approximated as a constant in [tn—1, tn] to im- 
493 prove the numerical stability in the case with a large time- 
494 Step. Though the linear approximation adopted in Eq. 33 
45 Shows higher accuracy to the constant approximation, the am- 
4» plitude is corrected with the solution of the EPKE. There- 
497 fore, the constant approximation can provide higher stability 
4» in large time-step cases without affecting the numerical accu- 
499 racy. In the PCQ-SCM solver, Eq. 34 is replaced by: 


491 


diis (tn) = Pvk.g (ca) gm (n) At, +wrvk,S,g(tn)Atn 


= dua (tn) 7 (A6 (56) 


500 


IH. C5G7-TD BENCHMARK PROBLEM 


501 


s2 The C5G7-TD transient benchmark is derived from the 
sos C5G7 benchmark model[31, 32]. It includes both 2D and 
so 3D models that are used to verify the transient capability for 
ss modeling a heterogeneous light water reactor without thermal 
sos feedback. In the following sections, the geometry, material 
so; and transient exercises in the CSG7-TD are described. 
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A. Geometry description 


The benchmark problem simulates a light water reactor 
with 8 uranium oxide (UO2) assemblies, 8 mixed oxide 
(MOX) assemblies, surrounded by the water reflector. The 
model is simplified into a quarter core using the reflective 
sts boundary on the north and west boundaries and the vacuum 
s514 boundary on the south and east boundaries, as shown in Fig. 
8-a. The four assemblies are numbered 1-4 respectively. The 
sie Size of the 2D model is 64.26 cm x 64.26 cm. The 3D model 
si; has the same planar layout as the 2D model, but in the axial 
direction, two identical water reflectors are added above and 
519 below the active region of the core as shown in Fig. 8-a. The 
s2 height of the 3D configuration is 171.36 cm while that of the 
active region is 128.52 cm. 


The UO2 assembly and MOX assembly have identical ge- 
ometry configurations, but different fuel pin compositions. 
Their sizes are 21.42 cm x 21.42 cm with a 17 x 17 pin 
layout. Each assembly consists of 289 pin cells, including 
264 fuel pins, 24 guide tubes, and 1 fission chamber as shown 
in Fig. 8-b. The UO2 assembly contains one type of fuel pin, 
while the MOX assembly contains three types of fuel pins 
with enrichment of 4.3%, 7.0% and 8.7%, respectively. 


e 
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511 
512 


515 


51 
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Each pin cell has a pitch of 1.26 cm and is simplified as 
two zones as shown in Fig. 8-c. Zone 1 is the circular area 
with the radius of 0.54 cm. It is filled with fuel, control rod 
or fission chamber. Zone 2 is located outside of the circular 
s3 zone and filled with moderator. 


536 B. Material description 


The C5G7-TD benchmark adopts multigroup macroscopic 

cross-sections and kinetics parameters (KPs) for eight materi- 
als: UO2 fuel, 7.0% MOX fuel, 4.3% MOX fuel, 8.7% MOX 
fuel, guide tube, fission chamber, moderator and control 
rod[32]. The cross-sections are provided in a 7-group struc- 
ture including transport-corrected total cross-sections Nig , 
absorption cross-sections Xag, fission cross-sections X fg, fis- 
sion spectra Xg, fission neutron yields v and scattering cross- 
sections 2£,,,;. The KPs are provide in an 8-delayed-group 
structure including neutron velocities vg, delayed neutron fis- 
sion spectra Xig, delayed neutron fractions 8; and decay con- 
stants A;. 
s For the 3D control rod movement problems, such as TD4 
sso problems, to reduce the rod cusping effects, the effective vol- 
ume fraction proposed by[6] is applied to homogenize the 
cross-sections of the guide tube and control rod in the par- 
tially rodded node. The effective volume fraction (y) of the 
ss4 rod is a polynomial function of the original volume fraction 
(x), which is the ratio of the inserted rod length to the node 
sss height. The polynomial function is given as: 


555 


y =0.3867848x + 0.1707878z? + 0.98878812? 
— 5.9535 775^ + 25.1805898a° — 63.184125249 
+ 98.38988022" — 92.49825962° + 48.0373368x? 
— 10.518205127° 


557 (57) 
sss Fig. 9 shows the relations between the original and effective 
volume fraction with Eq. 57 and linear assumption, respec- 


tively. 


559 


C. Transient description 


The C5G7-TD benchmark problem consists of six exer- 
cises labelled TDO-TD5. Each exercise includes a number 
of cases that simulate the movement of control rod banks or 
the change in moderator density. To model the control rod 
movements, the control rods in an assembly are referred to as 
a rod bank. Control rods in the same bank are simultaneously 
inserted or withdrawn. 

As illustrated in Table 1, TDO, TD1, and TD2 are 2D ex- 
ercises that require adjusting different control rod banks in 
a variety of ways. TDO contains 5 test problems. In TDO, 
control rods are abruptly inserted into the active region (10% 
of the core height) at the initial time and remain there for 1 
s7s S. Then, control rods are partially extracted (5% of the core 
sz height) and stay for 1 s. At the end of 2 s, the control rods are 
57 completely withdrawn from the active region. 
ss TDI and TD2 each contain 5 and 3 test problems. In 
s9 TD1 and TD2, control rods move linearly at a constant speed. 
sso From the starting point at t=0 s, control rods are inserted to 
ss: the maximum depth at 7-1 s and then returned to the start- 
see ing point at t2 s. The maximum inserted depth for TD1 and 
se» TD2 is 1% and 10% of the core height, respectively. In the 
se. 2D model, the control rod movement is simulated as the linear 
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Fig. 8. Geometry and composition of C5G7-TD benchmark 
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Fig. 9. Relation between original and effective rodded volume frac- 
tion 


substitution of the moderator-filled guide tube material with 
the control rod material in Zone 1 of the specified pin cells. 


TABLE 1. Scenarios of control rods movement defined in TDO, TDI 
and TD2 


Test problem TDO TD1 TD2 
1 Bank 1 Bank 1 Bank 1 
2 Bank 3 Bank 3 Bank 3 
3 Bank 4 Bank 4 N/A 
4 Bank 1, 3 and 4 Bank 1,3 and 4 N/A 
5 Bank 1-4 Bank 1-4 N/A 


TD3 is a 2D exercise simulating the reactivity insertion 
caused by the moderator density change in the reactor core. 
It has 4 test problems. It has 4 test problems. During the first 
second, the moderator density in all assemblies declines si- 
multaneously and at the same rate from the nominal value to 
the minimal value. After t=1 s, the density increases to the 
nominal value within the following second. These four test 
problems differ in the ratio of the minimum to the nominal 
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value as shown in 10. 
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Fig. 10. Fractional moderator density change in TD3 exercise 


TD4 is a 3D exercise that simulates the process of insertion 
and withdrawal of control rods. It contains 5 problems. In the 
initial phase, all control rods are placed in the top water re- 
flector, outside the active region. These five problems involve 
continual insertion and removal of various combinations of 
control rod banks. Fig. 11 depicts the control rod banks and 
the relative insertion length to the active core length for each 
problem. 


Similar to TD3, TD5 is a 3D transient event initiated by 
the moderator density change with all control rods fully with- 
drawal. The moderator density in the same assembly changes 
simultaneously. For each of the four test problems, the den- 
sity of the moderator varies in each of the four assemblies as 
shown in 12. 
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Fig. 11. Relative depth of control bank movement in TD4 exercise 
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IV. NUMERICAL RESULTS 


11 


617 the 2D problems and 48 processors for the 3D problems. 


Unless otherwise noted, all benchmark problems are sim- 
eie ulated with VITAS on a workstation with 16 processors for 
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Each fuel pin is meshed with three radial rings for the fuel 
zone, one radial ring for the moderator zone, and eight az- 
imuthal sectors, as determined by a preliminary h-p sensitiv- 
ity analysis[11]. In total, the mesh is comprised of 32 el- 
ements shown in Fig. 2. In the 3D calculation, the axial 
direction is discretized into 32 layers, and the thickness of 
each layer is 5.355 cm for TD5 exercise. In TD4 exercise, 
the axial direction is discretized into 48 layers to better de- 
scribe the control rod movement. For the spatial expansion 
within the node, the x-y plane is expanded with continuous 
FE basis functions of 2nd order, and the axial direction is ex- 
panded with orthonormal polynomials of 2nd order. In ad- 
dition, the interface is approximated with orthonormal poly- 
nomials of 2nd order. For the angular expansion, a 25x25 
square Legendre-Chebyshev cubature is used to evaluate an- 
gular integrals within the node, and the spherical harmonics 
functions of P23_P3 and P4 are implemented on the lateral 
and axial interfaces, respectively. Table 2 presents a sum- 
mary of the parameters. To reduce the computational effort, 
the 2D/1D approach (for 3D problem) and the quasi-reflected 
interface condition of up to P3 are applied. 


A. Steady-state results 


For steady-state calculations, Table 3 compares the eigen- ° 


value results of the 2D and 3D initial eigenvalue problems 
calculated by different solvers[5, 6]. VITAS results are con- 
sistent with those of other solvers. For the 2D case, the differ- 
ence in eigenvalues between VITAS and MCNPS is 29 pcm, 
which is less than the Monte Carlo uncertainty of 0.07%; for 
the 3D case, the difference between VITAS and PROTEUS- 
MOC is 25 pem. The minor deviations indicate that the model 
is accurate for the initial condition. In addition, the eigenval- 
ues of the forward and adjoint calculations are consistent with 
a difference of 1 pcm. 

Meanwhile, the normalized pin-wise fission rate distribu- 
tions of 2D and 3D problems are depicted in Fig. 13 and Fig. 


14, respectively. The normalization is implemented by scal- 
6 


ing the initial total fission rate to unity. 
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Fig. 13. Pin-wise fission rate distribution at the steady-state of 2D 
cases 
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Fig. 14. Pin-wise fission rate distribution at the steady-state of 3D 
cases 


B. Verification of direct-SCM solver 


To verify the direct-SCM solver in VITAS, VITAS re- 
sults are compared with results by MPACT[4] and PANDAS- 
MOC[6]. MPACT implements a novel Transient Multilevel 
(TML) approach that involves three levels of calculation: 3D- 
transport, 3D-Coarse Mesh Finite Difference (CMFD) and 
EPKE, in order to enhance computational efficiency. The 
coarse time-step solution of 3D-transport is corrected with 
fine time-step solutions of CMFD and EPKE. In PANDAS- 
MOC, the transient solution is solved by coupling the 2D 
MOC and 1D nodal methods accelerated by the multi-level 
CMFD methods. Due to the large reactivity insertion, the 
core power in the TD3-4 problem fluctuates rapidly. Conse- 
quently, the TD3-4 problem is chosen as a representative for 
the time-step size sensitivity analysis. 

The fractional core fission rates obtained with time-step 
sizes ranging from 10 ms to 200 ms are depicted in Fig. 15(a). 
Fig. 15(b) illustrates the reference solution and the percent- 
age differences between the reference solution and solutions 
with larger time-step sizes. The reference solution is obtained 
with a time-step size of 1 ms using the direct-SCM solver. 
The percentage difference e is calculated as: 


P (t) — P* (t) 


P(h) x 10096 


e(t) = (58) 
where the superscript * denotes the reference solution. Fig. 
15(a) indicates that as the time-step size increases to At — 
100 ms or 200 ms, there exists oscillation in the fission rate 
at the asymptotic stage (t = 2 ~ 3 s). Similarly, the os- 
cillation of the percentage difference is demonstrated in Fig. 
15(b). The numerical oscillation is attributed to the linear 
approximation of the amplitude frequency when a large time- 
step size is utilized. At the asymptotic stage, the control rods 
are completely withdrawn from the active region, and the fis- 
sion rate increases gradually due to the delayed neutrons from 
precursors. Due to the slower growth rate of the fission rate, 
the numerical oscillation becomes dominant at the asymptotic 
stage. It also explains that there is no observable oscillation 
when the fission rate changes rapidly during t = 0 ~ 2 s. 
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TABLE 2. Parameters of VITAS calculation model 


Model parameters 


Value 


Fuel pin mesh (fuel zone rings/moderator rings/azimuthal sectors) 3/1/8 


Volume spatial expansion order 

Surface spatial expansion order 

Volume angular integrals 

Surface Py order (lateral interfaces/axial interfaces) 
Number of axial meshes in 3D 


TABLE 3. Eigenvalue results for steady-states of 2D and 3D C5G7- 
TD cases 


Code Method 2D eigenvalue — 3D eigenvalue 
MCNP5 MC 1.18646+0.07% 
PROTEUS-MOC 3D MOC 1.18651 1.16469 
MPACT 2D/1D MOC 1.18666 1.16359 
NECP-X 2D/1D MOC 1.18695 

PANDAS-MOC 2D/1D MOC 1.18631 1.16512 
VITAS (Forward) VNM 1.18675 1.16494 
VITAS (Adjoint) VNM 1.18674 1.16493 


ess When the time-step size is decreased to 20 ms, the oscillation 
ex decays rapidly. In order to quantify the effect of the time- 
ess Step size further, Table 4 summarizes the highest percentage 
e»» differences with different time-step sizes. The RMS error is 
7 calculated using Eq. 59. According to Table 4, with a time- 
701 Step size of 20 ms, the solution agrees well with the refer- 
72 ence solution with an RMS error of less than 0.1%. Addition- 
703 ally, the benchmark suggests that the time-step size during the 
7 transient event should not exceed 25 ms during the transient 
75 event[32]. Therefore, 20 ms may be an effective time-step 
zæ Size for the direct-SCM solver to reduce the computational 
707 cost without losing the numerical resolution. If not otherwise 

, 708 Specified, the time-step size used by the direct-SCM solver in 
7» the following section is At = 20 ms. 


TABLE 4. Maximum difference and RMS error of different time- 
step sizes in TD3-4 


Step size (ms) Max . Diff. (70) RMS error (96) 


10 0.46396 0.035% 

20 0.826% 0.099% 

100 1.370% 0.499% 
a 200 2.889% 1.288% 
71 

N 2 
P (ty) — P* (t 

712 RMS = [ES | ( 2i ( ») (59) 
718 1. 2D results 


24 In this section, results of 2D exercises with ramp- 
71s perturbations are discussed, which include TD-1, TD-2 and 
ns TD-3. To access the accuracy of the simulation, the results of 
717 VITAS are compared to those of MPACT[4] and PANDAS- 
ns MOC[6]. Percentage differences are evaluated using MPACT 


719 
720 


p 


722 


72: 


© 


724 
725 
726 
727 
728 
729 
730 


73 


2 

2 

25x25 Square Legendre-Chebyshev cubature 
P. 23 -P. 3 / P. 3 

48 (TD4), 32 (TD5) 


1.0 


e 
to 


e 
o 


Relative fission rate 


e 
FS 


0.2 


0.0 0.5 1.0 1.5 2.0 2.5 3.0 
Time (s) 


(a) Fractional fission rate 


1.0 


0.8 


0.6 


o 
M 
Percent difference (96) 


Relative fission rate 


0.4 


M 
Y 


-- 10 ms-1ms 
-- 20ms-1ms 
-- 100 ms-1 ms 
== ims ---- 200 ms-1 ms 


0.2 


0.0 0.5 1.0 1.5 2.0 2.5 3.0 
Time (s) 


(b) Percentage difference 


Fig. 15. Fractional fission rate and percentage difference of differ- 
ence time-step sizes in TD3-4 


and PANDAS-MOC as references. The time-step sizes em- 
ployed in MPACT and PANDAS-MOC are 10 ms and 20 ms, 
respectively. It should also be noted that PANDAS-MOC uses 
an increased time-step size of 100 ms and 50 ms in the asymp- 
totic stage of TD1/2 and TD3 exercises. 

The fractional fission rates of TD1 and TD2 exercises are 
depicted in Fig. 16. In these problems, continuous changes of 
fractional fission rates are observed due to the movement of 
control rods with a constant speed. Despite that fact that the 
speed of insertion and withdrawal are identical, the increase 
rate of the fission rate after 1 s is less than the decrease rate 
before 1 s due to the reduction of the precursor concentra- 
tion. The fission rate then increases asymptotically to a value 


smaller than the initial value. Due to the deeper insertion of 
control rods in TD2, the fission rate changes more rapidly in 
TD2. Fig. 17 and Fig. 18 compare the direct-SCM solver's 
results with those of MPACT and PANDAS-MOC. Compar- 
isons reveal that VITAS with the direct-SCM agrees well with 
other codes for 2D control rods movement problems. Fig. 19 
depicts the fractional fission rate histories of TD3. Due to the 
similar pattern of perturbations caused by the moderator den- 
sity change, the fission rate exhibits similar tendencies. As 
seen in Fig. 20, the direct-SCM solver's results of TD3 are 
also consistent with those of other codes. 
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Fig. 16. Fractional fission rate of TD1 and TD2 exercises 


73 Table 5 summarized the RMS error and maximum percent- 
74 age differences in TD1, 2, and 3 exercises. Compared with 
the rest cases, the solutions of TD2-1 and TD3-4 show greater 
z deviations with maximum percentage differences of 1.0%. It 
747 is observed that, as depicted in Fig. 18 and Fig. 20, the max- 
imum differences appear in few time-steps after the introduc- 
tion of the reactivity. Since the reactivity insertion is substan- 
tial and fission rate changes rapidly in TD2-1 and TD3-4, dif- 
ferent temporal discretization schemes can lead to larger dis- 
752 crepancies. However, these discrepancies are less than 1.0%, 
7s3 Which implies that the direct-SCM solver of VITAS is capable 
754 Of performing accurate calculation for 2D ramp-perturbation 
755 problems. 
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Fig. 19. Fractional fission rate of TD3 exercise 


TABLE 5. Maximum difference and RMS error of 2D results 


Problem MPACT MPACT PANDAS PANDAS 
Max. Diff.(%) RMS Err.(%) Max. Diff.(%) RMS Err.(%) 
1-1 -0.13 0.08 0.14 0.06 
1-2 -0.16 0.08 -0.16 0.08 
1-3 -0.14 0.07 -0.14 0.06 
1-4 -0.17 0.07 0.17 0.06 
1-5 -0.17 0.08 0.27 0.08 
2-1 0.91 0.10 0.91 0.10 
2-2 -0.21 0.11 -0.23 0.15 
2-3 -0.12 0.05 -0.11 0.04 
3-1 0.25 0.11 0.29 0.16 
3-2 0.45 0.11 0.56 0.16 
3-3 0.72 0.11 0.72 0.15 
3-4 0.97 0.11 0.97 0.14 


2. 3D results 


756 


In this section, results of 3D TD4 and TDS exercises are 
discussed. In TD4 and TD5, VITAS employs time-step sizes 
of 25 ms and 20 ms, respectively. MPACT and PANDAS- 
MOC employ time-step sizes of 25 ms and 20 ms, respec- 
tively. At the asymptotic stage, the PANDAS-MOC time-step 
sizes are increased to 100 ms. 

The fractional fission rate history of TD4 exercise is de- 
picted in Fig. 21. In TD4-4 and TD4-5 problems, the fission 
765 rate is more complicated since multiple control rod banks are 
zæ inserted and withdrawn simultaneously. Fig. 22 compares 
767 VITAS results to those of MPACT and PANDAS-MOC. The 
z8 MPACT results are in well agreement with the other code re- 
zə sults, with VITAS results being the most comparable to the 
770 MPACT results. Oscillation in the percentage difference with 
71 à certain period can be observed in the stage of control rod 
72 movement. It is attributed to the different de-cusping tech- 
75 niques implemented in different codes. As shown in Table 
74 6, the maximum difference between MPACT and VITAS re- 
77 sults is -0.97%, whereas the maximum difference between 
77% PANDAS-MOC and VITAS results is 1.57%. The discrep- 
77 ancies between these codes in TD4 are caused not just by the 
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Fig. 17. Fractional fission rate comparisons of VITAS versus MPACT and VITAS versus PANDAS-MOC for TDI exercise 
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Fig. 18. Fractional fission rate comparisons of VITAS versus MPACT and VITAS versus PANDAS-MOC for TD2 exercise 


7: temporal nor spatial discretization scheme, but also by the de- 7e» linearly in different fuel assemblies, as shown in Fig. 


779 Cusping techniques. 


TABLE 6. Maximum difference and RMS error of TD4 
MPACT 


Problem 
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Max. Diff.(%) RMS Err.(%) Max. Diff.(%) RMS Err.(%) 


-0.97 
-0.77 
-0.77 
-0.75 
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MPACT 


0.28 
0.28 
0.29 
0.34 
0.39 


PANDAS 
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1.57 
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During the first two seconds, the fission rates decrease to the 
minimal values as the moderator densities decrease. In the 


785 Subsequent two seconds, the moderator densities are restored 
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Fig. 23 presents the evolution of the relative fission rates 


731 for the TD5 exercise, in which the moderator density varies 


to the nominal values, and the fission rate increases gradu- 
; ally. Fig. 24 compares the results of several codes. The 
comparison demonstrates that the direct-SCM solutions agree 
well with the MPACT and PANDAS-MOC results. As shown 
in Table 7, compared with the corresponding MPACT and 
PANDAS-MOC solutions in TD5-2 problem, the direct-SCM 
results show maximum differences of 0.87% and 0.85%, re- 
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Fig. 20. Fractional fission rate comparisons of VITAS versus MPACT and VITAS versus PANDAS-MOC for TD3 exercise 


795 3. Step-perturbation results 


1.0 


7e The TDO results of VITAS and MPACT are compared in 
797 Order to validate the direct-SCM solver for step-perturbation 
7 problems. Typically, the step-perturbation problem demands 
7» finer time-steps to capture the rapid change induced by the 
soo Step reactivity[7]. For a PCQM solver, a time-step is further 
discretized into several micro time-steps for EPKE solution. 
&2 The PCQM solver can therefore simulate the step reactivity 
s3 With a rather large time-step size without a significant loss of 
sos accuracy. However, the direct transient solver requires a finer 
à i à é 8 10 eos time-step size. To compare the solutions of the direct-SCM 

xo sos SOlver with those of MPACT in TDO exercise, variable time- 
so7 Steps are employed in VITAS. 


Fig. 21. Fractional fission rate of TD4 exercise sos In the TDO exercise, the reactivity is introduced abruptly at 
soo t = 0, 1,2 s, while the cross-sections stay constant for the re- 
eo maining transient activity. Consequently, a smaller time-step 

size can be employed immediately after the step reactivity 

si2 change; otherwise, a larger time-step size is used. The time- 

813 Step size employed in MPACT is 10 ms. Fig. 25 illustrates 

814 the variable time-step size employed in VITAS, which is 1 

sis ms during the first 20 ms following a step reactivity change, 
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TABLE 7. Maximum difference and RMS error of TD5 sis and 20 ms for the remainder of the transient process. 
Problem MPACT MPACT PANDAS PANDAS s7 The fractional core fission rates of TDO exercises are 
Max. Diff.(%) RMS Err.(%) Max. Diff.(%) RMS Err.(%)sıs shown in Fig. 26. The fission rates decrease at 0 s and 1 s as 
5-1 0.69 0.43 0.74 0.52 s19 the result of the abrupt insertion of control rods, then increase 
5-2 0.87 0.45 0.85 0.54 s2 promptly at 2 s with the withdrawal of control rods. After 
5-3 0.85 0.38 0.78 0.47 


821 2 S, all control rods are completely out of the active region, 


= pon po 0:34 0.12 22 and the fission rates hereby increase asymptotically to a value 
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Fig. 23. Fractional fission rate of TD5 exercise 


ses smaller than the initial value. Fig. 27 compares VITAS and 
s24 MPACT solutions, while Table 27 summarizes the error. The 
es maximum difference of all step-perturbation problems is ap- 
sæ proximately 0.21% in TDO-2 problem, while the differences 
827 in other problems are within 0.12%. They demonstrate that 
s28 the solutions of VITAS and MPACT are highly consistent, in- 
seg dicating that the direct-SCM solver is capable of generating 
a30 accurate solutions for step-perturbation problems. 
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Fig. 22. Fractional fission rate comparisons of VITAS versus MPACT and VITAS versus PANDAS-MOC for TD4 exercise 


TABLE 8. Maximum difference and RMS error of TDO 


Problem MPACT MPACT 
Max. Diff.(%) RMS Err.(%) 

0-1 -0.12 0.05 

0-2 -0.20 0.11 

0-3 -0.10 0.05 

0-4 -0.11 0.05 

0-5 -0.11 0.05 


Percent difference (%) 
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C. Performance of PCQ-SCM solver 


In the following section, the PCQ-SCM solver is tested 
against the direct-SCM solver using three problems: TD0-5, 
TD3-4, and TD5-1. The numerical accuracy and the compu- 
tational cost are applied as metrices for evaluating the perfor- 
mance of the method. 

Fig. 28 shows that the reactivity evaluated by VITAS is 
consistence with MPACT results. Fig. 29 compares the re- 
sults by direct-SCM and PCQ-SCM solvers for the TDO-5 
problem. The solution of direct-SCM solver is obtained with 
a constant time-step size of 1 ms, while the PCQ-SCM calcu- 
lations employ variable time-step sizes. To capture the step- 
reactivity, a step-size of 1 ms is employed in the first 1 ms 
following a step reactivity change, whereas 100 or 200 ms for 
the remaining. As shown in Fig. 29, the PCQ-SCM solu- 
tions are consistent with those of the direct-SCM solver, with 
the exception of the first 1 ms following the step reactivity 
change. As shown in Table 9, the percentage differences are 


co 
R 
© 


co 
a 
o 


eo 
E 


18 


VITAS 
— MPACT 
— PANDAS-MOC 


--- VITAS-MPACT TD5-1 
~~ VITAS-PANDAS TD5-1 |, 


VITAS 
—— MPACT 
—— PANDAS-MOC 


esset 


Relative fission rate 
S 
o 

Percent difference (%) 


i 
2 


--- VITAS-MPACT TD5-2 
== VITAS-PANDAS TDS-2 |. 


VITAS 
— MPACT 
—— PANDAS-MOC 


--- VITAS-MPACT TD5-3 
=> VITAS-PANDAS TD5-3 |. 


difference (94) 


Percent difference (26) 


Time (s) 


(a) TD5-1 


Time (s) 


(b) TD5-2 


2 3 
Time (s) 


(c) TD5-3 


VITAS 
—— MPACT 
—— PANDAS-MOC 


o 


Relative fission rate 


o 


--- VITAS-MPACT TD5-4 
== VITAS-PANDAS TDS-4 |. 


0.0 


b 
Percent difference (%) 


Time (s) 


(d) TD5-4 


Fig. 24. Fractional fission rate comparisons of VITAS versus MPACT and VITAS versus PANDAS-MOC for TD5 exercise 
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within 0.0896 in both time-step sizes. Note that when cal- 
culating the RMS error and maximum percentage difference 
displayed in Table 9, the results at the first 1 ms following the 
reactivity change are not included. The significant discrep- 
ancies following the reactivity change cause by the errors in 
the direct-SCM solution. The step size of 1 ms in the direct 
solver is insufficient to trace the step response of the fission 
rate; however, for the PCQ-SCM solver, the first 1 ms is fur- 
ther subdivided into 100 micro-steps. The micro-step size of 
0.01 ms is capable of simulating the step response accurately. 

In addition, Table 9 displays the computing time for differ- 
ent solvers. In comparison to the direct solution, the number 
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of TEVPs in the PCQ-SCM solution with a step size of 100 
ms decreases from 5000 to 53, and the computing time de- 
creases from 14 minutes to 7 minutes. The incommensurate 
decrease of computing time can be attributed to the fact that 
the computing time is dominated by response matrix forma- 
tion in 2D calculations. In the transient solvers of VITAS, the 
response matrices are only updated when the cross-sections 
are altered. In the TDO-5 problem, increasing time-step size 
to 100 ms accelerates the iteration for flux solving but not the 
response matrix formation. Therefore, the computational cost 
reduction brought by the PCQ-SCM solver is insignificant in 
TDO exercises. 
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Fig. 27. Fractional fission rate comparisons of VITAS versus MPACT and VITAS versus PANDAS-MOC for TDO exercise 


TABLE 9. Performance comparison of PCQ-SCM for TDO-5 problem (16 processors) 


Number of Max. Diff. RMS Err. Computing time 


Method At (ms) TEVPS (96) (96) (min) Speedup factor 
Direct 1 5000 N/A N/A 14 N/A 

PCQ 100 53 0.08 0.04 a 2.0 

PCQ 200 28 0.08 0.04 6 2.4 


Reactivity ($) 


Fig. 28. 
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Comparison of reactivity history for TD0-5 problem 


Fig. 30 illustrates the reactivity history of TD3-4 generated 


74 by VITAS and MPACT. Although the cross-sections vary lin- 
875 early in this problem, the reactivity curve versus time shows 
s76 non-linear. This behavior can be attributed to the spatial ef- 
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fects caused by the variation of the flux shape. Fig. 31 com- 
pares the solutions of direct-SCM and PCQ-SCM solver for 
the TD3-4 problem. The solutions of direct-SCM solver are 
obtained with fixed time-step sizes of 1 and 100 ms, and the 
PCQ-SCM calculations are performed with time-step sizes of 
100, 500 and 1000 ms, respectively. Fig. 31(b) shows that as 
the time-step size of the PCQ-SCM decreases, the PCQ-SCM 
solution gradually converges to the direct solution. The max- 
imum percentage difference in the solution with a step size 
of 100 ms is -0.23% as shown in Table 10. It is less than the 
direct solution with 10 ms shown in Table 4. Additionally, 
Fig. 31(b) also indicates that in the direct solution with 100 
ms, there are significant fluctuations at the asymptotic stage, 
which leads the maximum difference greater than 1.096. It 
indicates that the PCQ-SCM solver can accurately predict the 
fission rate and eliminate the numerical fluctuations with a 
large time-step size. 

In the meantime, the number of TEVP calculations for the 
PCQ-SCM calculation with a 100 ms step size is decreased 
from 10000 to 100, where the linear speedup is 100, and 
the computing time is reduced by a speedup factor of 78.9. 
The speedup of TEVP calculations is not proportional to the 
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Fig. 29. Comparison of direct-SCM and PCQ-SCM solutions for 
TDO-5 problem 


speedup of the total computing time. It is because that a 
larger time-step size makes it more difficult for amplitude fre- 
quency, fission source and flux more difficult to converge in 
the k — w iteration. However, the speedup of the overall com- 
puting time is close to the linear speedup of TEVP calcula- 
tions. It implies that the PCQ-SCM solver can reduce the 
computing time efficiently without a significant loss of the 
numerical accuracy in the TD3-4 problem. 

Fig. 32 shows the reactivity history of TD5-1 as evaluated 
by VITAS, MPACT and PANDAS-MOC. Note that the reac- 
tivity displayed in Fig. 32(a) is defined as p/ with the unit 
of $, whereas it is defined as p with the unit of pcm in Fig. 
32(b). Fig. 32(a) reveals that compared with the reactivity of 
MPACT and PANDAS-MOC, VITAS predicts an underesti- 
mated p/8. However, in Fig. 32(b), the reactivity p evaluated 
by VITAS is consistent with that of PANDAS-MOC. Con- 
sequently, the deviations in Fig. 32(a) can be attributed to 
the discrepancies of the effective delayed neutron fractions 3 
computed by different codes. 

Fig. 33 compares the results of direct-SCM and PCQ-SCM 
solvers for the TD5-1 problem. The direct solutions are ob- 
tained with time-step sizes of 1 and 100 ms. The macro 
time-step sizes used in the PCQ-SCM calculation are 100 
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Fig. 30. Comparison of reactivity history for TD3-4 problem 
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Fig. 31. Comparison of direct-SCM and PCQ-SCM solutions for 
TD3-4 problem 


and 500 ms. As shown in Fig. 33, during the first four sec- 
onds, the PCQ-SCM calculations provide an underestimated 
fission rate with respect to the direct-solution. When the 
cross-sections revert to the nominal value at the asymptotic 
stage, however, the fission rate values get closer to the ref- 
erence solution. The discrepancies are shown in Table 11. 
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TABLE 10. Performance comparison of PCQ-SCM for TD3-4 problem (8 processors) 


Number of Max. Diff. RMS Err. Computing time 


Method At (ms) TEVPS (96) (96) (min) Speedup factor 
Direct 1 10000 N/A N/A 6573 N/A 
PCQ 100 100 -0.23 0.12 83 78.9 
PCQ 500 20 -0.42 0.21 16 405.0 
PCQ 1000 10 -1.21 0.65 11 598.6 
0.0 Direct-SCM 1 ms 
1.0 * PCQ-SCM 100 ms 
* — PCQ-SCM 500 ms 
-0.5 *  Direct-SCM 100 ms 
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Fig. 32. Comparison of reactivity history for TD5-1 problem 


The maximum difference between the direct solution and the 
PCQ-SCM solution with 500 ms step size is -1.01%, and the 
RMS error is 0.66%. When the macro time-step size is re- 
fined to 100 ms, the solution and error show no significant 
change. Additionally, as shown in Fig. 33, compared with 
PCQ solutions, the direct solution with 100 ms shows smaller 
differences of less than 0.5%. The deviations between the 
direct- and PCQ-SCM solutions can be attributed to the in- 
consistent formulation of the kinetics parameters in VITAS 
transient models[5]. 


Meanwhile, Table 11 further summarizes the speedup of 
PCQ-SCM solutions. It shows that the speedup factors in 
TD5-1 problem are significantly lower than those for TD3- 
4 problem. The deviation in speedup is due to the different 


94: 


N 


943 
944 
945 
946 
947 
948 
949 
950 


95 


952 
953 
954 


955 


(a) Fractional fission rate 


Direct-SCM 1 ms -- PCQ error (100 ms) 
-- PCQ error (500 ms) 


-- Direct error (100 ms) 


= 
a 


Relative power 


9 
in 


Percent difference (96) 


e 
ES 


ud 
w 


Time (s) 


(b) Percentage difference 


Fig. 33. Comparison of direct-SCM and PCQ-SCM solutions for 
TD5-1 problem 


fractions of computing time in 2D and 3D problems. In VI- 
TAS, 99% of the overall runtime is consumed by obtaining 
response matrices and solving transport equation. Therefore, 
the following runtime analysis will focus on matrix formation 
and equation solving. Fig. 34(a) shows the fractional rumtime 
of matrix formation and equation solving for TD5-1 problem. 
When the time-step size is 1 ms, matrix formation and equa- 
tion solving account for 57.16% and 42.80% of the runtime, 
respectively. However, when using a time-step size of 100 
ms or 500 ms for PCQ-SCM solutions, equation solving is 
becoming dominant in the total runtime. Fig. 34(b) depicts 
the speedup of matrix formation, equation solving, and to- 
tal calculation. The speedup of the matrix formation is close 
to the linear speedup, whereas the speedup of the equation 
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TABLE 11. Performance comparison of PCQ-SCM for TD5-1 problem (48 processors) 


Number of Max. Diff. RMS Err. Computing time 


Method At (ms) TEVPS (96) (96) (min) Speedup factor 
Direct 1 500 N/A N/A 6820 N/A 
PCQ 100 50 -1.00 0.66 271 25.2 
PCQ 500 10 -1.01 0.65 99 68.9 
solving is significantly below the linear speedup. As noted oes lation of the direct-SCM, the SCM transforms the transient 


previously, the fission source and flux require more iterations 
to get converged when time-step size is increased. Therefore, 
the equation solving speedup is not proportional to the linear 
speedup. As the time-step size increases, the runtime time is 
dominated by the equation solving, and the speedup of the to- 
tal calculation is primarily determined by the speedup of the 
equation solving. 
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transport equation into a TEVP, and then the solution of the 
TEVP is solved by the heterogeneous VNM. In addition, the 
PCQ-SCM is proposed and developed to reduce the compu- 
tational cost. The predicted neutron flux in the PCQ-SCM is 
obtained by solving the TVEP, and is then corrected by the 
EPKE solution. For the formulation of the EPKE, the meth- 
ods for computing the reactivity and adjoint flux within the 
heterogeneous VNM are discussed. Due to the second-order 
even-parity form of the transport equation in the VNM, the 
adjoint equation can be easily solved by creating the adjoint 
source without generating additional response matrices. For 
the calculation of reactivity, the exact reactivity is determined 
using the dynamic eigenvalue acquired from the TEVP, with- 
out the need to reconstruct the angular flux or apply perturba- 
tion theory. 

To verify the direct-SCM solver, four C5G7-TD 2D bench- 
mark exercises and two 3D benchmark exercises are exam- 
ined. The direct-SCM solver's transient results are in good 
agreement with those previously published by MPACT and 
PANDAS-MOC. In 2D problems, the percentage difference 
relative to other reference solutions is within 0.9796, while 
in 3D problems, the difference is within 0.97% and 1.57% 
compared with MPACT and PANDAS-MCC results, respec- 
tively. It demonstrates the correctness of the direct-SCM 
implemented in VITAS. In addition, the performance of the 
PCQ-SCM solver is compared to the performance of the 
direct-SCM solution on TD0-5, TD3-4, and TD5-1 problems. 
Due to the response matrix update strategy implemented in 
VITAS, the PCQ-SCM can minimize the computing time for 
step-perturbation problem TDO-5 by approximately twofold 
as it decreases the number of iterations with larger time steps. 
For ramp-perturbation problems, PCQ-SCM with larger time- 
steps is capable of drastically lowering the number of re- 
sponse matrices and power iterations. The computing time 
for various problems, such as TD3-4, can be saved by one to 
two magnitude without sacrificing accuracy. Meanwhile, the 
analysis to the computing time reveals that the speedup effi- 
ciency of PCQ-SCM is determined by the runtime fractions 
of the matrix formation and equation solving. In 2D ramp- 
perturbation problem, where the matrix formation typically 
dominates the runtime, the PCQ-SCM has a high speedup ef- 
ficiency; whereas in 3D problems, it exhibits a lower speedup 
efficiency since the equation solving is dominant in the run- 
time. 

Future efforts will focus on improving the parallel ef- 
ficiency of the code so that it can be deployed on clus- 


This work discusses the direct-SCM and PCQ-SCM tran- 105 ters to solve more complex problems. Moreover, the high- 
ee sient models of the VITAS code, as well as the verification of 1v6 fidelity transient multi-physics coupling of the neutron trans- 
s; these models using the C5G7-TD benchmark. In the formu- 1017 port method with thermal-hydraulic simulations, and the un- 
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1018 Certainty analysis[33] to the transient multi-physics calcula- 1023 appeared to influence the work reported in this paper. 
1019 tion will be investigated. 
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